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The path integral molecular dynamics (PIMD) method provides a convenient way to compute the 
quantum mechanical structural and thermodynamic properties of condensed phase systems at the 
expense of introducing an additional set of high-frequency normal modes on top of the physical vi- 
brations of the system. Efficiently sampling such a wide range of frequencies provides a considerable 
thermostatting challenge. Here we introduce a simple stochastic path integral Langevin equation 
(PILE) thermostat which exploits an analytic knowledge of the free path integral normal mode 
frequencies. We also apply a recently-developed colored-noise thermostat based on a generalized 
Langevin equation (GLE), which automatically achieves a similar, frequency-optimized sampling. 
The sampling efficiencies of these thermostats are compared with that of the more conventional 
Nose-Hoover chain (NHC) thermostat for a number of physically relevant properties of the liquid 
water and hydrogen-in-palladium systems. In nearly every case, the new PILE thermostat is found 
to perform just as well as the NHC thermostat while allowing for a computationally more efficient 
implementation. The GLE thermostat also proves to be very robust delivering a near-optimum sam- 
pling efficiency in all of the cases considered. We suspect that these simple stochastic thermostats 
will therefore find useful application in many future PIMD simulations. 



I. INTRODUCTION 

Atomistic computer simulations are now routinely used 
to shed light on the behavior of condensed phase systems 
of interest in physics, chemistry, biology and materials 
science. Typically these simulations assume that the nu- 
clei can be treated as classical particles, even when the 
interactions between them are obtained using ah initio 
methods. For systems comprised of heavy atoms at high 
temperatures, this is often quite a reasonable assump- 
tion. However, whenever any light atoms such as hydro- 
gen are present at room temperature or below, quantum 
mechanical zero point energy and tunneling effects in the 
nuclear motion can play an important role in determining 
the properties of the system. 

Since the 1980s it has been possible to include quantum 
mechanical effects in simulations of structural and ther- 
modynamic properties such as radial distribution func- 
tions and average potential and kin etic energies by per- 
forming imaginary time path integral^ (pj) simulations. 
These simulations exploit the isomorphism between the 
quantum mechanical partition function of the physical 
system and the classical partition function of an extended 
problem consisting of a necklace of replicas of the system 
connected by harmonic springs P 

Early implementations of this approach used Monte 
Carlo (MC) methods to sample the configuration space, 4 
but these methods often proved to be inefficient unless 
carefully constructed collective moves were used. It was 
then suggested by Parrinello and RahmarJ^l that mo- 
menta could be assigned to the replicas, allowing the 
sampling to be performed by molecular dynamics (MD). 
In principle, this allows many of the techniques originally 



developed for use in classical MD simulations to be used 
in a path integral context. 

In practice, the stiff harmonic springs between the 
replicas lead to inefficient and nonergodic dynamics when 
microcanonical trajectories are used to generate ensem- 
ble averages.^ A simple stochastic Andersen thermostat- 
ting scheme^ was originally suggested to overcome this 
problemP However, stochastic approaches have since 
largely been abandoned following the development of de- 
terministic Nose-Hoover chain thermostatsp^n^l which 
have now become the gold standard for performing 
path integral molecular dynamics (PIMD) simulations.^ 
These thermostats generate ergodic, canonical averages 
and provide a conserved quantity which can be used to 
check the integration time step, at the expense of intro- 
ducing sets of auxiliary chain varia bles which add to the 
complexity of the calculation! 12 * 13 ! 

In the last few years stochastic approaches have re- 
gained attention in a variety of different contexts, and 
many of their drawbacks have been overcome. For exam- 
ple, a strategy to define a conserved quantity has been 
developed which not only allows one to check the qual- 
ity of the integration but also allows for the correction 
of sampling errors due to the use of a finite time steppe 
Indeed the total energy of the system minus the accu- 
mulated heat absorbed from the thermostat will clearly 
be a conserved quantity for any thermostatting scheme. 
Moreover new stochastic methods have been int roduc ed 
which show very promising sampling properties! 15 * 16 ! In 
particular, a stochastic velocity rescaling thermostat has 
been developed which couples to the total kinetic en- 
ergy of the system rather than the individual momenta 
of the particles within itP^ This allows canonical sam- 
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pling to be achieved with a much smaller disturbance 
to the Hamiltonian trajectory. In addition, it has been 
shown how a generalized (colored-noise) Langevin equa- 
tion thermostat can be tuned to sample a ver y wid e range 
of frequencies simultaneously and efficiently! 16 * 17 ! 

The purpose of the present paper is to exploit these 
recent developments by introducing two new stochastic 
path integral thermostats that are competitive in terms 
of sampling efficiency with the Nose-Hoover chain ther- 
mostat but simpler to implement and cheaper to use. 
To demonstrate this, we have performed an extensive 
set of benchmark calculations, computing the statistical 
sampling efficiencies of various thermostats for a number 
of physical observables in two very different condensed 
phase systems: a flexible model of liquid water at room 
temperature and the diffusion of an interstitial hydrogen 
atom in palladium. 

The outline of the paper is as follows. Section II be- 
gins by reviewing the theory of PIMD and discussing how 
the equations of motion for the path integral replicas (or 
ring polymer beads) can be integrated in the absence 
of a thermostat. This sets the scene for the remain- 
der of the section, which introduces two new stochas- 
tic thermostats for path integral simulations, compares 
them with the well-established Nose-Hoover chain ther- 
mostat, and ends with a brief discussion about the dif- 
ference between global and local thermostatting. Section 
III presents the results of our calculations on the liquid 
water and hydrogen-in-palladium systems and Section IV 
presents our conclusions. 



II. THEORY 

A. Path integral molecular dynamics 

Consider a system of distinguishable particles de- 
scribed by a Cartesian Hamiltonian of the form 

N 2 

H = J2^ + V (^--->^ (!) 

i— 1 1 

in which the potential energy V(q\, . . . , q^) is such that 
the quantum mechanical partition function 

Z = ix [e~^ H ] (2) 

is well defi ned (with j3 — 1/k-oT). After a standard Trot- 
ter produclP^l] discretization of the trace, this partition 
function can be written as 

where / = Nn and f3 n = f3/n. Here H n (p, q) is the clas- 
sical Hamiltonian of a fictitious ring polymer consisting of 



n copies of the system connected by harmonic springsp^ 
J ff„(p,q)=^(p,q) + V4(q), (4) 

where 

N " (\ 0')l2 i \ 

i=l j=l \ 1 / 

(5) 

and 

n 

V n (d) = ^2V(q[ j \...,q^), (6) 

with w„ = 1/ j3 n h and qf = qf . The error in Eq. (3) is 
0(1/ n 2 ) and so vanishes in the limit as n —t oo. 

The path integral molecular dynamics (PIMD) 
methooP uses this classical isomorphism as a computa- 
tional tool to calculate quantum mechanical equilibrium 
properties of the form 

(A) = Itr [e-^A] . (7) 

For example, the potential energy of the system is given 
by 

{v) - j^Wz I dfp I d/ q e ^" ff " (p ' q)v »(q)> ( g ) 

where the estimator V n (q) involves an average over the 
beads of the ring polymer necklace, 

V»(q) = -K(q) = V(q[ 3 \ . . . , q%\ (9) 

n n _. =i 

and the kinetic energy can be calculated as 

where 7^(q) is most efficiently chosen to be the centroid 
virial estimator] 20 ! 21 ! 

(11) 

with 

1 " 

3=1 

The errors in Eqs. (8) and (10) are again 0(l/n 2 ) and 
may therefore be neglected for sufficiently large n. 

It is clear from these equations that the potential and 
kinetic energies are obtained from estimators which do 
not involve any reference to the ring polymer momenta, 



3 



and indeed one can easily integrate out the momenta 
from Eqs. (8) and (10) to leave purely configurational 
averages. The momen ta w ere originally introduced by 
Parrinello and Rahman-* as a sampling tool to facili- 
tate the exploration of configuration space by molecular 
dynamics. As we have written it, the free ring polymer 
Hamiltonian H®(p, q) in Eq. (5) corresponds to one par- 
ticular choice of the Parrinello-Rahman mass matrix in 
which the physical particle mass is assigned to each ring 
polymer bead. This is the choice that is used in the ring 
polymer molecular dynamics (RPMD) appro ximation to 
real time quantum correlation functions ,1221221 although of 
course real-time properties are immediately lost as soon 
as any thermostat is switched on, as we shall do in this 
paper. 

The PIMD method is based on the observation that, 
since the classical ring polymer trajectory 



Pt 



dH n (p t ,q t ) 
dq t 



and q t 



dH n (p t ,q t ) 
dp t 



(13) 



conserves both the Boltzmann factor e ^H n {p t ,it) = 
e -^„ff„(po,qo) anc j £ ne phase space volume element 
d f p t d f q t = d f p a dfq Q , Eqs. (8) and (10) can be re- 
written as 



(A) 



1 



(2TTh)fZ 



d f p I d f 



A n (r t ). 



This implies that static equilibrium properties such as the 
potential and kinetic energies can be obtained by time- 
averaging along ring polymer trajectories whose initial 
conditions are sampled from the Boltzmann distribution 



/o(Po,qo) 



1 



(2irh)fZ 



-/3„i?„(po,qo) 



(15) 



Alternatively, and more efficiently, one can combine the 
sampling with the time evolution by attaching a thermo- 
stat to the ring polymer dynamics as we shall describe 
below. Before we do this, however, let us first explain 
how to integrate the equations of motion in Eq. (13) in 
the absence of a thermostat. 



where L = Lq + Ly is the Liouvillian associated with 
H n (p, q) and Lo and Ly are those associated with 
#°(p,q) andK(q). 

The exact evolution generated by i?"(p, q) is simplified 
by transforming the ring polymer from the bead repre- 
sentation to the normal mode representation, 



# } = £*< Wc k and fff > = 5>(«C, fcl 



(17) 



3=1 



3=1 



where in the case of even n the elements of the orthogonal 
transformation matrix C are 



Cj k - 



fl/n, k = 

/2/n cos (2njk/n) , 1 < k < n/2 - 1 
^ljn(-l) j , k = n/2 

^2~Jn sin (2Trjk/n) , n/2 + l<k<n-l. 

(18) 

In the normal mode representation, H^(p, q) becomes 



fl!(p,q)=EE' fe 



i=l k=0 



2mi 



;miUJ k [q, 



2r^( fe )i2 



with 



LOk = 2u n sin (kir/n) . 



(19) 



(20) 



In view of this, the algorithm implied by Eq. (16) for 
integrating the equations of motion in Eq. (13) through 
a time interval At is as follows: 



At dV{q\ 



(J) 



Oq 



(J) 



(21) 



if^E^- f ^E«- ( 22 ) 

3=1 3=1 



cosu>kAt — rriiUJk sinajfc A<\ Ip 



J ■<- y[l/miU) k ]smu k At cosw fe Ai J J ' 

(23) 



B. Ring polymer time evolution 

A convenient way to integrate Eq. (13) is based on the 
splitting of the Hamiltonian in Eq. (4) into a sum of a 
free ring polymer part H®(p, q) and a potential energy 
part V n (q). This suggests combining the exact evolutions 
generated by these two parts of the Hamiltonian in a 
symplectic integration scheme in which the phase space 
density evolves under the influence of the symmetric split 
operator propagator 



-AtL 



e -(At/2)i v e -AtL 0e -(At/2)L v 



(16) 



ri— 1 — l 

ft^EC*^. ?f ^£<^#\ (24) 

k=0 k=0 



AtdV{q[ 3 \. 



Oq 



(J) 



(25) 



The first step in this algorithm is an exact evolution of 
the ring polymer momenta through a time interval At/2 
under the influence of the Hamiltonian V n (q). The sec- 
ond is a transformation from the bead representation to 
the normal mode representation. The third is an exact 
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evolution of the ring polymer coordinates and momenta 
through a time interval At under the influence of the 
free ring polymer Hamiltonian H®(p,q). The fourth is 
a transformation from the normal mode representation 
back to the bead representation, and the fifth is a fur- 
ther evolution of the ring polymer momenta through a 
time interval At/2 under the influence of the Hamilto- 
nian V n (q). 

Because the algorithm consists of a sequence of exact 
evolutions under the influen ce of a pproximate Hamiltoni- 
ans, it is exactly sy mplect ic 1211251 w hich implies that the 
ring polymer phase space volume will be conserved for 
any time step At. However, the algorithm will not be 
very accurate (for example, the ring polymer Hamilto- 
nian will not be very well conserved) unless At is suffi- 
ciently small. When there is only one ring polymer bead, 
the transformations in Eqs. (22) and (24) become redun- 
dant, and Eqs. (21), (23) and (25) reduce to the standard 
(second order) velocity Verlet method 2 ^ for integrating a 
classical trajectory. It is clear from Eq. (16) that the algo- 
rithm remains accurate to 0(At 2 ) for any number of ring 
polymer beads. Note also that, in view of the nature of 
the orthogonal transformation matrix C in Eq. (18), the 
transformations in Eqs. (22) and (24) can be done using 
fast Fourier transform routines. These are so efficient 
that we have not bothered to minimize the number of 
transformations between the bead and normal mode rep- 
resentations in the implementation of the algorithm we 
have given above and shall also not bother to do so in the 
implementation of the thermostats described below|23 



C. A path integral Langevin equation thermostat 

Bussi and Parrinello have recently explained how a 
simple (white noise) Langevin thermostat can be com- 
bined with the velocity Verlet algorithm to give an effi- 
cient sampling of th e can onical distribution in classical 
statistical mechanics ! 14 * 28 ! Since PIMD is simply classi- 
cal molecular dynamics in an extended phase space (al- 
beit with a canonical distribution at n times the physical 
temperature), and since the algorithm we have given in 
Eqs. (21) to (25) is a direct generalization of the veloc- 
ity Verlet algorithm, it is straightforward to adapt this 
Langevin thermostat to the present context and use it to 
thermostat PIMD. 

This can be done in two distinct ways, depending 
on whether one chooses to thermostat the ring poly- 
mer beads or the normal modes. Choosing the latter for 
reasons that will become apparent below, one replaces 
Eq. (16) with 



-AtL 



e -(At/2)i T g-(Af/2)iv e -AtLo e -(At/2)L v g -(At/2)i 7 _ 



(26) 

where L 7 is the part of the Liouvillian L = Lq + Ly + L y 
in the Fokker-Planck equation for the Langevin phase 
space density that introduces the friction and the ther- 
mal noise Bussi and Parrinello have shown that this is 



equivalent to adding the following algorithmic steps both 
before and after Eqs. (21) to (25) P 



Pi ^ A^Pi °J fc ' 



~(fc) . (fc) ~(k) . 

p\ <- c\ 'p\ + 




n-1 
fe=0 



> C jk pf\ 



(27) 



(28) 



(29) 



Here £^ is an independent Gaussian number (a normal 
deviate with zero mean and unit variance) that is differ- 
ent for each physical degree of freedom, each ring polymer 
normal mode, and each invocation of Eq. (28), and the 
coefficients and (4 are^l 



„(*) 



„(*0 



-(At/2)7 (fc) 



(fc)l 



= V 1 -[4' 



(30) 



(31) 



All that remains to complete the specification of this 
thermostat is to specify the normal mode friction coef- 
ficients "f^ k \ The advantage of working in the normal 
mode representation is that these friction coefficients can 
be chosen to give an optimal sampling of the canonical 
distribution for the free ring polymer, in a sense that we 
shall now explain. 

In the normal mode representation, the Langevin dy- 
namics of each mode of the free ring polymer is that of 
an uncoupled harmonic oscillator, 



dt % 



~(k) 



dt. 



~(fc) 
Pi 



2 ~(fe) 



7 (fc)p(*) 



(32) 

where Q (t) represents an uncorrelated, Gaussian- 
distributed random force with unit variance and zero 
mean = and (^ fc) (0)^ (fc) (t)) = S(t).} In view 

of this, the optimum friction coefficient 7W will be that 
which gives the smallest autocorrelation time of the har- 
monic oscillator Hamiltonian and thus the most rapid 
(Boltzmann-wcightcd) exploration of free ring polymer 
energy shells. The relevant autocorrelation time 

th = — 1 — 2 £ ((H(0) (H))(H(t) (H))) dt 

(33) 



5 



where 



^( fc )l2 



H 



2 771, 



:miUJ k [q 



(34) 



can be worked ou t anal ytically for the Langevin dynamics 
in Eq. (32) and isPlMI 



th = 



1 

7*) 



7 (fe) 



(35) 



The optimum value of 7W for the excited (fc > 0) modes 
of the ring polymer is therefore = 2cjfc. However, this 
prescription cannot be used to thermostat the centroid 
mode, as it gives 7^ = 2wo = (purely Hamiltonian dy- 
namics without a thermostat). We shall therefore define 
a separate thermostat time constant To for the centroid 
mode and specify the normal mode friction coefficients 
in Eq. (30) as follows: 



v( fe ) = 



1/t , k = 
2u> kl k > 0. 



(36) 



With this specification, the algorithm in Eqs. (27) 
to (29) provides a very simple way to thermostat a PIMD 
simulation that only requires a single input parameter 
(to) and that has been tuned to sample the internal 
modes of the ring polymer as efficiently as possible. The 
simplicity arises because the tuning is based on the fre- 
quencies of the free ring polymer internal modes, which 
are known analytically and are independent of the inter- 
actions in the system (and therefore transferrable from 
one system to another). The price to be paid for this 
is that the optimum friction coefficients for the free ring 
polymer may not be quite the same as those for the inter- 
acting ring polymer, although we would expect them to 
be very similar for the highest frequency internal modes 
which are well separated from the vibrations of the phys- 
ical system. Just how well this "path integral Langevin 
equation" (PILE) thermostat works in practice will be 
investigated for two different systems in Sec. III. 



wide range of frequencies, including both the physical vi- 
brations of the system and the high frequency ring poly- 
mer internal modes. That it is possible to construct such 
a GLE thermostat is shown in Fig. 1, which compares the 
sampling efficiency of a colored noise thermostat with 
that of a white noise thermostat with a friction coeffi- 
cient of 7 = luq. The sampling efficiency is defined^ 
as k(uj) = [ojtv , where Ty(uS) is the autocorrela- 
tion time for the potential energy of a harmonic oscil- 
lator with frequency uj. This provides an indication of 
how efficiently the thermostat explores the thermally ac- 
cessible configuration space of a vibration with this fre- 
quency. The white noise thermostat is optimally efficient 
(k = 1) when u) = ujq, but its efficiency decreases lin- 
early for higher and lower values of w such that an effi- 
ciency of k > 0.2 is only obtained for a frequency range 
spanning two orders of magnitude. By contrast, the col- 
ored noise thermostat achieves k > 0.2 for a frequency 
range spanning more than four orders of magnitude. In 
a 32-bead path integral simulation of liquid water at 300 
K, this should be enough to sample everything from the 
highest frequency internal mode of the ring polymer at 
2n/f3hc ~ 13,300 cm -1 down to the diffusive modes of 
the liquid at frequencies as low as 1 cm . 

The colored noise thermostat in Fig. 1 is one of a family 
of thermostats that have recently been constructed^ by 
exploiting the equivalence^ between the non-Markovian 
dynamics of the GLE and Markovian dynamics in a 
higher dimensional space. These thermostats have a sim- 
ple matrix form that is straightforward to implement on 
top of the ring polymer evolution algorithm in Eqs. (21) 
to (25). Choosing to work in the bead representation 
rather than the normal mode representation, since there 
is no advantage in this case in making the transformation 
to normal modes, one simply replaces Eqs. (27) to (29) 
witrP 1 



(37) 



Here £^ is a vector of n s + 1 independent Gaussian num- 
bers, 



D. A generalized Langevin equation thermostat 

Another recent development in thermostatting has 
been the construction of colored-noise Langevin ther- 
mostats for variety of different problems, ranging from 
separately thermostatting the electrons and the ions in 
Car-Parrinello molecular dynamics^ to efficiently gen- 
crating configurations consistent with the quantum me- 
chanical canonical ensemble without introducing any 
path integral beadsPH The key point here is that the 
generalized Langevin equation (GLE) has sufficient flex- 
ibility in its memory kernel and its colored noise to be 
optimized for a variety of different applications. 

In the present PIMD context, we would like to choose 
the colored noise so as to efficiently thermostat a very 




(38) 



is a vector containing the momentum of bead j and n s 
auxilliary momenta , and 



Ci = e 



-(At/2)-, T 



(39) 



and 



QfC 2 =I-Cfd, (40) 



are the appropriate matrix generalisations of Eqs. (30) 
and (31). 

Given a suitable friction matrix 7, the implementation 
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of this algorithm is as follows. One begins by construct- 
ing Ci from Eq. (39) and C 2 from Eq. (40). Ci can 
be constructed by combining a low-order Taylor series 



expansion of e ( At / 2 



with P matrix-squaring op- 



erations, and C2 can be constructed by Cholesky decom- 
position of I — CfCi. These matrices are then stored 
for use in each iteration of Eq. (37), which requires two 
matrix- vector multiplications and the generation of n s + 1 
normal deviates for each degree of freedom i and each ring 
polymer bead j. 

In most large-scale applications, the effort required by 
this algorithm will be negligible compared with the eval- 
uation of the forces in Eqs. (21) and (25). If it is not, 
one can exploit the fact that the canonical distribution 
is invariant under the action of Eq. (37) for any time 
step At, and only apply the thermostat once in every m 
time steps with At in Eq. (39) replaced by m At !^ With 
this modification, the algorithm clearly provides a com- 
putationally very tractable way to thermostat a PIMD 
simulation. 

All that remains is to specify the friction matrix 7 in 
Eq. (39). In the calculations described below, we used 
n« = 4 and 



7 = ui A, 



(41) 



where the elements of the (5 x 5) matrix A are given in 
Table I. These matrix elements were obtained in an au- 
tomated procedure^ designed to optimize the efficiency 
k(u>) of the thermostat, which we have already compared 
with that of the corresponding white noise Langevin ther- 
mostat (n s = 0, 7 = ujq) in Fig. 1. In that comparison, 
loq was chosen arbitrarily, but in actual applications it 
is a physical input parameter that can be used to tune 
the response of the thermostat to a particular frequency 
range, which one sees to be roughly 0.01 w < w < 100 oj 
from the results in Fig. 1. One can of course equally well 
replace this input parameter with a thermostat time con- 
stant To = l/2tJo, which we shall do in our comparison 
of various thermostats in Sec. III. 

For the sake of comparison, we will also perform some 
calculations in Sec. Ill with a white noise Langevin equa- 
tion thermostat with friction 7 = uj q = l/2r . This white 
noise Langevin scheme is closely related to an Andersen 
thermostat in which the momenta are resampled with a 
probability jAt after each time step At. The two meth- 
ods have a similar sampling efficiency in all the examples 
we have studied, and for this reason the results with the 
Andersen thermostat will not be reported. 



E. The Nose-Hoover chain thermostat 

As discussed in the Introduction, the gold standard 
against which to compare the above stochastic PILE and 
GLE thermostats is a deterministic Nose-Hoover chairP^ 
applied separately to each physical degree of freedom and 
each ring polymer bead (or ring polymer normal mode) 



TABLE I: Dimensionless elements Aij of the GLE friction 
matrix A in Eq. (41). This matrix was constructed using 
the options n a = 4, range = 10 4 , and "fit for potential 
energy", using the freely available software on the website 
|http:/ / gle4md.berli os.de/] 



At 



1 1 2.468046483820e+l 

1 2 3.618484148135e-2 

1 3 1.529754837748e+0 

1 4 -4.832976901522c+0 

1 5 3.075592122514c+l 

2 1 -3.690906142217e-2 
2 2 1.140757569304e-5 
2 3 9.580998002948e-2 
2 4 -2.633785831010e-2 

2 5 5.628596350432e-2 

3 1 -1.967695128248e+0 
3 2 -9.580998002948e-2 
3 3 1.803797247061e-l 
3 4 6.834981703810e-l 

3 5 -1.326536043516e+0 

4 1 -1.376606646573e+0 
4 2 2.633785831010e-2 
4 3 -6.834981703810e-l 
4 4 3.538593762043e+0 

4 5 1.527314768745c+0 

5 1 2.893495089306c+l 
5 2 -5.628596350432e-2 
5 3 1.326536043516e+0 
5 4 -1.527314768745e+0 
5 5 4.108827095695c+l 



One convenient way to do this that reduces to a 
recommended^! operator splitting in the classical limit 
is to replace Eq. (16) with 

e -AtL ^ e -(At/2)L NH c e -(At/2)Lv e -AtL 0e -{At/2)L v e -(At/2)L NHC 

where Lj^hc is the part of the Liouvillian L — L + 
Ly + Lnhc for the extended system dynamics that in- 
volves the Nose-Hoover chain variables.^ Choosing to 
work in the normal mode representation as in the case 
of the PILE thermostat, the net effect of this operator 
splitting is that Eq. (28) is replaced by the evolution of 
the following system of non-linear differential equations^ 
through a time interval of At/ 2: 



d k ) - 



dt. 



_(fe) 
-Pi 



n {k) 



(43) 
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10 10 2 10 3 10 4 10 5 

cj [cm -1 ] 



FIG. 1: The vibrational spectrum from a path integral molecular dynamics simulation of a flexible water model (blue). In order 
to include the internal modes of the ring polymer, the density of states has been computed from the velocity autocorrelation 
function of the individual ring polymer beads. This vibrational spectrum is compared with the harmonic-limit sampling 
efficiency Ky of a white-noise Langevin thermostat (dashed) and an optimal GLE thermostat (red). The latter has been 
optimized so as to give constant sampling efficiency over a frequency range spanning four orders of magnitude, centered 
geometrically on ujq = 1/2tq. The diffusion coefficient for the water model, and that computed for the GLE dynamics in the 
free-particle limit, are also shown in terms of the "diffusion frequency" u>r> = k^T/mD. 



dt ' { m, p n ^QW 



-.(*) 

(k) n i,2 



(44) 



d ik) 



_(*) 

4? 5^ ( 4 5) 



Q(fe) 



kilt] 2 1 \ 



dt 



CO 



dt 



(46) 



(47) 



Here n^f and rj^f are the momentum and position vari- 
ables of the Nose-Hoover chain attached to the fc-th nor- 
mal mode of the ring polymer in the i-th degree of free- 
dom, for I = 1, ... , L. 



The system of differential equations in Eqs. (43) to (47) 
has been discussed extensively in the literature on Nose- 
Hoover chains and numerous numerica l methods have 
been proposed for integrating it! 13 * 32 * 33 ] The accuracy of 
the integration through each time interval At/2 can be 
checked separately for each i and k by monitoring the 
locally conserved quantity 
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and the accuracy of the overall symmetric split operator 
propagator in Eq. (42) can be checked by monitoring the 



globally conserved quantity 



jv n-i l 
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In order to compare this thermostat with the stochas- 
tic thermostats introduced above, we need to specify the 
normal mode thermostat masses Q {k \ These play an 
analogous role to the normal mode friction coefficients 
7^) in the PILE thermostat introduced in Sec. II. C. 
The Nose-Hoover chain masses that give the most effi- 
cient sampling of the canonical distribution for the free 
ring polymer are^ QV°) = l/p n uj\. However, we cannot 
use this prescription to thermostat the centroid mode, as 
luq = 0. As in the case of the PILE thermostat, we can get 
around this difficulty by defining a separate thermostat 
time constant tq for the centroid mode and replacing ojq 
with I /2to , so that the normal mode masses that appear 
in Eqs. (43) to (47) are specified as follows: 



Q 



(k) _ 
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fci 



k = 
k > 0. 



(50) 



This completes the specification of the NHC thermo- 
stat that we shall use to compare with the stochastic 
PILE and GLE thermostats in Sec. III. Several other 
groups have proposed slightly different implementations 
of the NHC thermostat for PIMD simulations, but we 
would not expect the differences to have any major im- 
pact on our results. For example, the original implemen- 
tation by Tuckerman et a?P^ involved transforming to 
"staging" variables rather than normal mode variables, 
and employed a different operator splitting than the one 
we have used in Eq. (42). Staging variables are a valid al- 
ternative to normal mode variables and there are a ram- 



ber of different ways in which one can do the operator 
splitting, all of which should in principle converge on the 
same result. We have used normal mode variables here to 
make the connection with the free ring polymer evolution 
algorithm in Eqs. (21) to (25), and chosen the operator 
splitting in Eq. (42) to emphasize the connection with the 
PILE thermostat in Eqs. (27) to (29). A general feature 
of NHC thermostats for path integrals is that they require 
the solution of a large number of systems of non-linear 
differential equations of the form in Eqs. (43) to (47). 
Since these differential equations must be solved numer- 
ically these thermostats are more complicated than the 
stochastic thermostats we have described above. 



and 



sign [a\ = sign 
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(53) 



(54) 



and c = e~ 7 . We shall call this global version of the 
PILE thermostat PILE-G and the local version in which 
the centroid is thermostatted using Eq. (28) PILE-L in 
the numerical comparisons presented below. 



F. "Global" versus "local" thermostatting 

All of the thermostats we have discussed so far are "lo- 
cal", in the sense that each degree of freedom and each 
ring polymer bead (or normal mode) is separately ther- 
mostatted. This is the simplest way to ensure that the 
kinetic and potential energies of every single particle in 
the system are thermalized as rapidly as possible, and it 
is therefore likely to be the most efficient way to sam- 
ple local properties such as the energy of an interstitial 
hydrogen atom in palladium. On the other hand, there 
are many properties that are sensitive to the slow, col- 
lective motions of all of the atoms in the system, such as 
the dielectric constant of liquid water. Local thermostats 
will almost certainly be too aggressive to give an efficient 
sampling of these collective properties because with even 
a moderate amount of friction they tend to to inhibit 
the Hamiltonian diffusion. There has therefore been a 
renewed interest recently in the development of "global" 
thermostats that are attached to the system in a more 
gentle way and have a less disruptive effect on the diffu- 
sion. 

In particular, Bussi et a/! 15 * 34 ! have derived a global 
version of the finite time step Langevin propagator in 
Eq. (28) that acts as a thermostat on the kinetic energy 
of all N degrees of freedom rather than on each degree 
of freedom separately. When implementing this global 
thermostat in a path integral context, it is desirable to 
apply it only to the centroid mode, and to leave more 
aggressive local thermostats attached to the less ergodic 7 
excited ring polymer internal modes. In practice, this 
amounts to replacing just the centroid (k = 0) component 
of Eq. (28) with the following stochastic velocity rescaling 
algorithmPMl 



-(0) 1 ~(o) 
Pi <~ a Pi i 



(51) 



where 



a2 = :l (i-c)([err+Et 2 [cH z ) , 2C (o) . i c(i - C ) 



2/3 n K 



2(3 n K ' 
(52) 



The Nose-Hoover chain thermostat can also be applied 
globally and indeed this is how classical canonical en- 
semble simulations are often carried out. Choosing again 
to keep the thermostatting of the excited ring polymer 
modes local, the modification amounts in this case to 
replacing the centroid (k = 0) components of Eqs. (43) 
to (47) witfP 1 



dt 



~(0) 
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dp 
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(55) 



(56) 



(57) 



(58) 



(59) 



and to making the appropriate modifications to Eqs. (48) 
and (49); for example Eq. (48) is replaced bjB^ 



N 

E 

=i 



2m, 
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(60) 



Notice in particular that the N separate Nose-Hoover 
chains in the local thermostat have been replaced with a 
single chain that is coupled to the total kinetic energy of 
the ring polymer centroid via Eq. (56). We shall call this 
global version of the NHC thermostat NHC-G and the 
local version in which the centroid is thermostatted using 
Eqs. (43) to (47) NHC-L in the comparisons presented 
below. 
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One could clearly also imagine trying to develop a 
global version of the GLE thermostat introduced in 
Sec. II. D. However, little benefit would be expected from 
this generalization, because the capability of the GLE to 
adapt automatically to the vibrational modes present in 
the system relies on the fact that an independent thermo- 
stat is applied to each degree of freedom. In any event, 
the two global thermostats we have just described (PILE- 
G and NHC-G) are sufficient for us to illustrate when it 
is preferable to use a global thermostat rather than a lo- 
cal thermostat for path integral simulations, as we shall 
do in the following section. 



III. RESULTS AND DISCUSSION 

It is well known that molecular dynamics can be used 
to efficiently explore regions of phase space where the free 
energy surface contains local minima separated by bar- 
riers comparable to the average thermal energy. Ther- 
mostats can further enhance this sampling, particularly 
in the case of PIMD where poorly ergodic, harmonic 
necklace modes are present. In some cases inefficient 
sampling can also arise when a statistically significant 
portion of phase space consists of regions separated by 
high free energy barriers. In these situations, which are 
not the concern of the present work, the fact that PIMD 
is simply classical molecular dynamics in an extended 
phase space should in principle enable one to adopt many 
of the accelerated dynamics techniques that have been 
developed to speed up the rates of transitions across bar- 
riers in classical simulations. 

When this issue is not present a quantitative way assess 
the sampling efficiency is to consider the correlation time 
of an observable A, 

ta = — 1 — 2 £ ((A(0) (A))(A(t) (A))) dt, 

(61) 

in which the angular brackets denote an average of the 
appropriate path integral estimator _4„(q) along a ther- 
mostatted PIMD trajectory. This is related to the sta- 
tistical uncertainty ca in the expectation value of the 
observable (A) by 



eaocW- — , (62) 

V ^sim 

where t slm is the total simulation time. Therefore, one 
would like to make ta as small as possible so as to reduce 
the uncertainty in (A) for a simulation of a given length. 
In the previous section we have used the correlation time 
of the total energy for a harmonic oscillator as a tentative 
measure of the efficiency of sampling for a Langevin equa- 
tion [see Eqs. (33) and (35)]. Here we will instead per- 
form PIMD simulations of two typical condensed-phase 
systems and compute the correlation times of a number of 
different physical observables so as to compare the sam- 



pling efficiencies of various thermostats in more realistic 
(anharmonic) situations. 



A. Liquid water 

Liquid water is ubiquitous in biological systems and 
important in many applications in chemistry and ma- 
terials science. Nuclear quantum effects have a signifi- 
cant influence on its properties and it has therefore been 
a common target of PIMD studies ever since the first 
path integral s imul ations of the liquid were performed in 
the mid 1980sP^U ft \ s a \ so a highly structured liquid 
in which changes in the local structure depend on com- 
plex collective rearrangements of the hydrogen bonding 
network, which are difficult to sample. Because of this, 
liquid water is a highly relevant test case for our compar- 
ison. 

Simulations were performed using the q-TIP4P /F flex- 
ible, four-site water potential which has recently been pa- 
rameterized to reproduce many of the static and dynami- 
cal properties of water in path integral simulations.^ For 
each thermostat and value of the time constant r , we ran 
a trajectory for 12 ns at a temperature of 298 K and a 
density of 0.997 g cm~ 3 . The simulation box contained 
216 molecules, with each atom described by a 32 bead 
ring polymer. The equations of motion were integrated 
as described in Sec. II, a nd the ring-polymer contraction 
scheme of Markland et a/! 38 * 39 ! was used to accelerate the 
evaluation of the long-range electrostatic interactions. 

In order to assess the efficiencies of the different ther- 
mostats, we first computed the correlation time tt of the 
centroid virial estimator of the kinetic energy. While for 
a system of classical particles the kinetic energy is simply 
distributed according to the Boltzmann law, for quantum 
systems there is also a position-dependent contribution 
which requires one to sample the forces experienced by 
the internal modes of the ring polymer. This contribution 
is included in the centroid virial estimator in Eq. (11), 
which can be written equivalently in the normal mode 
representation as 

rn(q) = ^-^Ei;>-^(q), (63) 
H i=l k=l 

where Pj; k \q) = —dV n (q)/dq^ is the force on the k-th 
excited (k > 0) normal mode of the ring polymer that 
arises from the physical interaction potential. 

Fig. 2 compares the kinetic energy correlation times tt 
for a white noise Langevin equation (WNLE) themostat, 
the GLE themostat, and global and local versions of the 
PILE and NHC thermostats as a function of the ther- 
mostat time constant To. The performance of the simple 
WNLE thermostat, which yields efficient thermalization 
in the high friction limit but is very inefficient for large 
To, shows that effective sampling of 7^(q) requires strong 
coupling to the high frequency necklace modes. This 
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FIG. 2: Correlation times for the centroid virial kinetic en- 
ergy (left) and the potential energy (right) obtained from path 
integral molecular dynamics simulations of liquid water. The 
six panels show the dependence of these quantities on the 
thermostat relaxation time to for the WNLE and GLE ther- 
mostats and for local and global versions of the PILE and 
NHC thermostats. 



FIG. 3: Correlation time for the squared dipole moment of 
the simulation box (left) and the computed molecular centre- 
of-mass diffusion coefficient (right) obtained from path inte- 
gral simulations of liquid water. The six panels show the 
dependence of these quantities on the thermostat relaxation 
time to for the WNLE and GLE thermostats and for local 
and global versions of the PILE and NHC thermostats. 



is achieved automatically in the PILE and NHC ther- 
mostats, which specifically target the ring polymer nor- 
mal modes, overcome the ergodicity problems observed in 
microcanonical PIMD and consistently yield a low value 
of tt- Note that this is true regardless of the method 
(global or local) that is used to thermostat the centroid, 
which does not contribute to the kinetic energy estima- 
tor in Eq. (63). The quality of the GLE thermostat also 
deteriorates more slowly than that of the WNLE ther- 
mostat with increasing To, and is constant as long as all 
of the ring polymer necklace modes are included in the 
optimal range of fitted frequencies (see Fig. 1). 

For comparison with tt, we also report in Fig. 2 the 
correlation time of the potential energy, Ty. This is much 
more sensitive to the softer modes, and to the diffusive 
motion of the centroid in particular. Efficiently sam- 
pling the potential energy is considerably (at least an 
order of magnitude) more difficult than sampling the ki- 
netic energy, as can be seen from the correlation times 
in the figure. This is because one must carefully bal- 
ance the strong coupling needed to ensure ergodicity of 
the polymers with the need to avoid overdamping which 
causes the thermostat to become an additional bottle- 



neck to diffusion. The diffusion bottleneck is quite ap- 
parent in Fig. 2 from the significant increase in Ty which 
is observed for all of the local thermostats when t is 
decreased below around 100 fs. 

A reasonable efficiency in the sampling of the poten- 
tial energy can be achieved by carefully tuning the sim- 
ple WNLE thermostat, with the optimal coupling being 
To ~ 500 fs. However, tuning the thermostat in this way 
for every system one might wish to simulate is clearly un- 
desirable and large increases in the potential energy auto- 
correlation time result from a less than optimal choice of 
the WNLE time constant. In contrast, the GLE, which 
also requires no a priori knowledge of the internal fre- 
quencies of the ring polymer, obtains a nearly constant 
sampling efficiency (ry ~ 2 ps) over the five orders of 
magnitude of To we have considered. 

Turning now to the targeted schemes which exploit a 
knowledge of the ring polymer normal mode frequencies, 
one sees from Fig. 2 that the potential energy is sam- 
pled most efficiently by PILE-G and NHC-G, which use a 
global thermostat for the centroid. These schemes partic- 
ularly excel in the case of strong coupling (i.e., for a small 
relaxation time tq), implying that a rapid rescaling of 
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the total (centroid) kinetic energy helps to speed up the 
canonical sampling of the potential energy. At the same 
time, since the global thermostats are coupled to the to- 
tal kinetic energy rather than to individual momenta, the 
trajectory and the dynamical properties are only slightly 
disturbed,^ and hence the ability of Hamiltonian dy- 
namics to generate complex, collective rearrangements 
can be fully exploited. This is clear from the superiority 
of the global thermostats over the local thermostats for 
this problem. An important final point from Fig. 2 is 
that the simple stochastic PILE-G thermostat provides 
equally rapid sampling of both the kinetic energy and the 
potential energy of liquid water as the NHC-G thermo- 
stat, without the need to evolve any extended variables. 

The superiority of global thermostats over local ther- 
mostats for liquid water is further supported by the re- 
sults in Fig. 3, where we show the correlation time of the 
squared dipole moment of the supercell and the center-of- 
mass diffusion coefficient obtained from the thermostat- 
ted trajectories. The squared dipole moment is related to 
the dielectric constant and is known to converge slowly 
as it requires the reorientation of many individual water 
molecules, which in turn requires a concerted rearrange- 
ment of the hydrogen-bonding network. In this case, 
microcanonical dynamics is extremely effective, and the 
best sampling is obtained with very weak coupling to the 
thermostats. For the local WNLE, PILE-L and NHC- 
L thermostats, there is a clear correspondence in Fig. 3 
between the reduction of the diffusion coefficient due to 
over-damped dynamics and the degradation of T42 as the 
thermostat time constant tq is decreased. The global 
PILE-G and NHC-G thermostats, and to a lesser extent 
also the GLE thermostat, give much more effective sam- 
pling of the squared dipole moment for small values of 
To and also have a less disruptive effect on the diffusion. 
The properties of the correlated noise we have used in 
the GLE are clearly such that its effect on diffusion is 
much less severe than that of the WNLE, 1 ^ despite its 
strong coupling to the high frequency necklace modes. 
As a consequence, its effects are intermediate between 
those of local and global schemes. 

Combining the results in Figs. 2 and 3, we see that the 
PILE-G and NHC-G thermostats are the most efficient 
at sampling the properties of liquid water that we have 
considered, with the GLE also doing rather well (and cer- 
tainly much better than the simple WNLE). However, 
one cannot conclude from this that it will always be bet- 
ter to use a global rather than a local thermostatting 
scheme. Local schemes enforce a canonical distribution 
on each of the individual components of the momentum, 
and therefore ensure a homogeneous distribution of en- 
ergy throughout the system. Global schemes only mon- 
itor the overall temperature of the system, and rely on 
anharmonic coupling between different regions and dif- 
ferent frequency ranges to reach equipartition. For these 
reasons, one should be particularly careful when using 
global thermostats, as the inefficient sampling of partic- 
ular internal degrees of freedom may be masked by the 



very efficient equilibration of the total temperature. A 
few local properties should therefore also be monitored to 
ensure that equilibrium has been reached, and in general 
a local thermostat is to be preferred whenever inhomoge- 
neous or quasi-harmonic problems are treated, as we will 
demonstrate in the next subsection. 



B. Hydrogen in palladium 

The second example system we shall consider is both 
inhomogeneous and quasi-harmonic: the motion of an 
atomic hydrogen (H) interstitial in the lattice of palla- 
dium (Pd). In this more harmonic system, strong anti- 
correlations are present in the underdamped limit. In 
some cases this can be beneficial, and methods have 
been devised which explicitly enhance anti-correlations.^ 
However, our aim here is to estimate how efficiently 
uncorrelated configurations are generated by the ther- 
mostatted dynamics. An analysis based on correlation 
times of the form in Eq. (61) would therefore be mislead- 
ing, as anti-correlations can significantly reduce ta and 
mask the presence of very long relaxation times. For this 
reason, the correlation times reported in this subsection 
are computed as the integral of the absolute value of the 
normalized autocorrelation function, 

ta = — 1 — a f |<(A(0) - (A))(A(t) - (A)))\ alt. 
{A ) - [A) Jo 

(64) 

The calculations used a supercell containing 256 Pd 
atoms and a single H atom. A 10 bead ring polymer 
was used for H while, in view of their large mass, the Pd 
atoms were treated as classical particles. The Pd lattice 
parameter was 3.89 A and the temperature was chosen to 
be 350 K, where diffusion of the H atom is still relatively 
fast (D = 0.16 A 2 ps _1 ). Interactions were modeled by 
an embedded atom potential.^ We performed 8 ns of 
simulation with a time step of 0.5 fs for each choice of 
thermostat and To- 

The thermostatting of the classical metal lattice and 
the H ring polymer were performed separately. The same 
coupling parameters and thermostatting method were 
used for the metal and for the centroid of the H atom 
necklace. With this setup global thermostats on the H 
atom centroid act on only three degrees of freedom. Nev- 
ertheless, significant differences are still observed with 
respect to fully local schemes. 

In Fig. 4 we plot the correlation times of the radius 
of gyration of the ring polymer and of the H atom ki- 
netic energy. The path integral estimators for both of 
these quantities depend solely on the internal necklace 
modes, and the methods which target these modes specif- 
ically (PILE and NHC) offer clear advantages. In par- 
ticular, the radius of gyration of the H atom ring poly- 
mer, which is almost completely decoupled from the slow 
motion of the ring polymer centroid, is seen to be ther- 
mostatted equally well by PILE and NHC. This clearly 
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FIG. 4: Integrals of the absolute values of the normalized au- 
tocorrelation functions f (Eq. 64) for the ring polymer radius 
of gyration (left) and the kinetic energy (right) of a hydrogen 
atom interstitial in a palladium lattice, obtained from path in- 
tegral molecular dynamics simulations. The six panels show 
the dependence of these quantities on the thermostat relax- 
ation time to for the WNLE and GLE thermostats and for 
local and global versions of the PILE and NHC thermostats. 



FIG. 5: Integral of the absolute value of the normalized au- 
tocorrelation function f (Eq. 64) for the potential energy of 
a hydrogen atom interstitial in a palladium lattice (left) and 
the diffusion coefficient of the hydrogen atom (right) obtained 
from path integral molecular dynamics simulations. The six 
panels show the dependence of these quantities on the thermo- 
stat relaxation time to for the WNLE and GLE thermostats 
and for local and global versions of the PILE and NHC ther- 
mostats. 



demonstrates that, by targeting the normal modes with 
either thermostat (stochastic or deterministic), the er- 
godicity problem of PIMD is completely resolved. The 
GLE is also seen to give rapid sampling of the H atom 
radius of gyration and kinetic energy without exploiting 
any knowledge of the normal mode vibrational frequen- 
cies, unless its fitted range is shifted so that it does not 
encompass all of the relevant spectrum (tq > 100 fs). 

Let us now consider two quantities which are more sen- 
sitive to the physical forces on the H atom: its contribu- 
tion to the potential energy (defined as the energy of the 
H in Pd system minus the energy of the Pd lattice at the 
same configuration) and its diffusion coefficient. The cor- 
relation time Ty of the H atom potential energy and the 
diffusion coefficient D obtained from the thermostatted 
PIMD simulations are shown in Fig. 5, where the dif- 
ferences between this problem and liquid water become 
more apparent. Global schemes are now very inefficient, 
as the frequency mismatch between the H atom motion 
and the Pd lattice motion results in very slow equilibra- 
tion unless a separate thermostat is coupled to each H 
atom centroid degree of freedom. This is clear from the 



results of the PILE-G and NHC-G thermostats, which 
do not perform any better for any value of the relaxation 
time t than the corresponding local thermostats perform 
in the limit of weak coupling (tq = 10 4 fs). We expect 
that this problem would be exacerbated even further if 
several interstitial H atoms were present and a global 
thermostat were to be attached to the total kinetic en- 
ergy of their centroids, as we did in applying the global 
schemes to liquid water. 

For the local schemes Ty responds sharply to changes 
in the thermostat relaxation time tq. This occurs because 
the range of H atom vibrational frequencies involved is 
relatively narrow, and by choosing an appropriate fric- 
tion it is possible to obtain nearly-optimal sampling even 
with a simple WNLE. The GLE thermostat is almost in- 
dependent of the choice of t , at the expense of a slight 
decrease in efficiency with respect to the optimal white- 
noise friction. This behavior is consistent with the ana- 
lytical estimates of the efficiencies of the two thermostats 
in the harmonic limit (Fig. 1), which show the GLE to 
have a lower maximum efficiency in exchange for the far 
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broader range of frequencies for which it gives near opti- 
mal sampling. 

In contrast to the water calculations we see that the 
diffusion coefficient of the H interstitial in Pd is less af- 
fected by overdamping, which manifests itself in Fig. 5 
only when a very small value of To is used. Unlike water, 
where diffusion is sensitive to orientational modes which 
lie at higher frequencies in the spectrum, the diffusion 
mechanism of H atoms in Pd is a much simpler lattice 
mediated process. The dynamics can therefore tolerate 
a higher level of disturbance from the thermostat before 
being hampered. 

Note finally from Fig. 5 that the deterministic NHC-L 
thermostat does not suffer from such a dramatic degra- 
dation in sampling efficiency for small r as its stochas- 
tic counter-part PILE-L. Indeed, it has previously been 
observed^ that for purely harmonic potentials the NHC 
dynamics yields nearly-optimal efficiency for all frequen- 
cies smaller than uq = 1/ \J j3 n = 1/2tq. Since the 
NHC equations lack rotational invariance, however, their 
performance depends critically on the relative orienta- 
tions of the eigenvectors of the Hessian with respect to 
the directions to which the thermostat is applied.^ The 
optimal behavior is only obtained with the correct align- 
ment and the sampling efficiency degrades significantly 
for other alignments in an anisotropic potential, as we 
have observed here in the case of liquid water (compare 
the NHC-L results for fy in Fig. 5 with those for ry in 
Fig. 2 in the limit of small r ). 

IV. CONCLUDING REMARKS 

Deterministic NHC thermostats have become the de 
facto standard for the difficult problem of canonical sam- 
pling in path integral molecular dynamics. Stochastic 
thermostats are however considerably simpler, and they 
provide a physically very appealing way to model the in- 
teraction of a system with a heat bath. Recent advances 
in the application of stochastic methods to various molec- 
ular dynamics sampling problems have therefore encour- 
aged us to reconsider these methods as an alternative to 
Nose-Hoover chains for path integral simulations, as we 
have done in this paper. 

Overall, we believe that the results we have obtained 
are very promising. For two distinctly different physical 
systems - room temperature liquid water and a hydro- 
gen atom interstitial in a palladium lattice - the simple 
stochastic PILE thermostat that we have introduced per- 
forms just as well as the deterministic NHC thermostat 
for nearly every property we have considered (see Figs 
2-5). The only exception we have found is for the po- 



tential energy of a H atom in a Pd lattice in the strong 
coupling limit (small To), where the sampling efficiency of 
the NHC thermostat degrades more gently (see the bot- 
tom left hand panel in Fig. 5). In all of the other cases we 
have considered there is essentially no difference in terms 
of sampling efficiency between our stochastic PILE and 
a deterministic NHC. 

The GLE thermostat we have investigated has also 
been found to perform very well, especially given that 
(unlike the NHC and PILE thermostats) it does not ex- 
ploit an analytic knowledge of the free ring polymer nor- 
mal modes. A particularly nice feature of this thermo- 
stat is that it combines a sufficiently strong coupling to 
the internal modes of the ring polymer to ensure ergodic 
dynamics with a more gentle perturbation on the low 
frequency centroid motion that only mildly inhibits the 
Hamiltonian diffusion. As a result of this feature, the 
GLE behaves in many respects as a compromise between 
a global and a local thermostat, and provided its time 
constant To is chosen such that the range of fitted fre- 
quencies encompasses the entire spectral range of interest 
it gives close to optimal sampling in every situation!^ 

One final comment we should make when assessing 
the relative merits of these various thermostats concerns 
their computational cost. The implementation of the 
Langevin equation is clearly very simple in both its white 
noise (scalar) and colored noise (matrix) forms, because 
an exact propagator can be obtained for the linear free- 
particle LE. In contrast, a multiple time step scheme 
is mandatory when solving the non-linear Nose-Hoover 
chain equations if one wants to avoid a drift in the con- 
served quantity. In an ab initio PIMD simulation, where 
the computational effort is dominated by the evaluation 
of the physical forces acting on each ring polymer bead, 
the extra effort that is required to solve the Nose-Hoover 
chain equations will often be negligible. However, in 
the present simulations the evaluation of the forces was 
made comparatively inexpensive by the use of empirical 
force fi elds a nd a highly efficient ring polymer contraction 
scheme J2S129J ^ s a consequence, the different thermostats 
were found to have quite a significant impact on the com- 
putational cost of the calculations.^ 

In summary, we firmly believe that stochastic methods 
should be reconsidered for use in path integral molecular 
dynamics (and MD in general) , as their sampling efficien- 
cies are comparable to that of the most commonly used 
deterministic scheme and they have a number of practical 
advantages. We are certainly now using Langevin equa- 
tion thermostats in our own PIMD simulations and we 
expect that the results we have presented in this paper 
will encourage others to do so also. 
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